AO-A032  046 
UNCLASSIFIED 


NAVAL  RESEARCH  LAB  WASHINGTON  D C F/G  20/14 

NUMERICAL  SOLUTION  OF  THE  PARABOLIC  EQUATION  WITH  THE  PEP  CODE.(U) 

OCT  76  B E MCDONALD 

NRL-MR-3385  NL 


\<f\ 

ADA03P046 

m 

END 

DATE 

FILMED 

1 -77 

ADA0  32046 


Numerical  Solution  of  the 
Parabolic  Equation  with  the  PEP  Code 


B.  E.  McDonald 

Plasma  Dynamics  Branch 
Plasma  Physics  Division 


October  1976 


NAVAL  RESEARCH  LABORATORY 
Washington,  D.C. 


Approved  for  public  release;  distribution  unlimited 


SECURITY  CLASSIFICATION  of  This  PAGE  '*han  Data  Entarad) 


||Hrr^ 

I EES; 


REPORT  DOCUMENTATION  PAGE 


NUMERICAL  SOLUTION  OF  THE  PARABOLIC 
EQUATION  WITH  THE  PEP  CODE  e 


L 1 

j E.jAlcDonald 


I PERFORMING  ORGAN!  IATiON  NAME  ANO  ADDRESS 


Nava)  Research  Laboratory 
Washington,  D.C.  20375 


II.  CONTROLLING  OFFICE  NAME  ANO  AOORESS 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


J.  TYPE  OF  REPOST  A PERIOD  COVEREO 

Interim  report  on  a continuing 
NRL  problem. 


• performing  org,  report  number 


• CONTRACT  or  grant  numbered 


is.  security  class.  (•*  tma  raport) 

UNCLASSIFIED 


is«.  OECLASSl  FI  CATION/  DOWNGRADING 

schedule 


i«  distribution  statemcnt  rof  mil  Rapon; 


Approved  for  public  release;  distribution  unlimited.  ^ h I ^ ^ ^ v>  ^ 


17  DISTRIBUTION  STATEMENT  (of  f/»*  abatract  antarad  In  Block  70,  II  dlllarant  from  Report) 


it.  Supplementary  notes 


IS  KEY  WOROS  (Conffnu#  on  ravaraa  alda  ll  nacaaaary  and  Idantlty  by  block  numbar) 


Propagation 
Random  media 
Refraction 
Parabolic  equation 


20  ABSTRACT  (Contlnua  on  ravaraa  alda  II  nacaaaary  and  Idantlty  by  block  numbar) 

We-prepose  a centered  leapfrog  difference  scheme  for  solution  of  the  parabolic  equation  for  ( 
wave  propagation  in  random  media.  We  stability  and  resolution  requirements  for  the  scheme.  * 
Results  of  two  numerical  integrations  are  in  good  agreement  with  those  obtained  by  others  using 
different  methods. 

Cn. 


DD  , 1473  COITION  OF  I NOV  «9  It  OBSOLETE  j 

S/N  0I02-014-  ««01 


SECURITY  classification  OF  this  PAGE  rRTl.il  Daia 


$ 61  °i50 


CONTENTS 

1.  THE  WAVE  EQUATION  1 

2.  THE  PARABOLIC  APPROXIMATION  3 

3.  NUMERICAL  INTEGRATION  OF  THE  PARABOLIC 

EQUATION  4 

A.  FINITE  DIFFERENCE  SCHEME  4 

B.  STABILITY  ANALYSIS 5 

C.  BOUNDARY  CONDITIONS  8 

D.  RESULTS  FROM  THE  PEP  CODE 9 

SUMMARY 11 

REFERENCES  13 


iii 


f 


NUMERICAL  SOLUTION  OF  THE  PARABOLIC  EQUATION 
WITH  THE  PEP  CODE 

1.  THE  WAVE  EQUATION 

We  are  concerned  with  obtaining  numerical  solutions  to  the  three- 
dimensional  wave  equation 

[V2  + k 2 (1  + y(r))]  E(r)  = 0 (1 ) 

in  a region  of  space  where  the  index  of  refraction  varies  slowly  with 
r.  Subject  to  certain  assumptions,  (1 ) describes  radio  propagation 
through  a spatially  varying  plasma*,  or  acoustic  wave  propagation 
through  a medium  in  which  the  sound  speed  is  a function  of  space2  . 

The  solution  E is  assumed  to  have  a time  dependence  eht , and  the  wave- 
number  kQ  is  gu/c,  where  c is  the  speed  of  light  for  the  radio 
application,  or  a reference  sound  speed  for  the  acoustic  application. 

For  radio  wave  propagation,  Ei'r)  is  a scalar  component  of  the 
electric  field,  and  y(r)  is  the  known  electric  susceptibility  of  the 
plasma.  Equation  (l)  neglects  a depolarization  term  whose  magnitude 
is  approximately  kQ  |Evy|  for  small  y.  If  X varies  on  a length  scale 
L,  the  depolarization  term  is  smaller  than  the  kQ2  y term  of  (1 ) by  a 
factor  ~/kQL.  We  shall  be  interested  in  applications  in  which 
L » ^/k  , so  the  depolarization  term  is  expected  to  be  negligible. 

The  use  of  a scalar  y also  assumes  implicitly  that  magnetization  is 
unimportant;  namely,  that  <u  is  much  greater  than  the  electron  cyclo- 
tron frequency. 

For  acoustic  applications,  E is  the  sound  pressure,  and  y is  a 
measure  of  the  sound  speed  fluctuation  such  that  the  local  sound  speed 
Note:  Manuscript  submitted  September  20,  1976. 


1 


IS  c 


hk.  + v(r)  . Just  as  in  the  radio  wave  application,  a term  involv- 


ing the  gradient  of  v has  been  neglected,  but  again,  its  magnitude  is 

down  from  that  of  the  yk  2 term  by  a factor  ^/k  L . 

* o o 

For  solutions  of  (l),  let  us  make  the  substitution 


E(r ) = u(r) 


(2) 

(3) 


ij"  kfz  )dz 

&i.r;  = u j e , 

where 

k(z)  = kQ  /l  + XQ  z ) , 
and  y (z)  is  some  "background"  x profile  along  an  assumed  propagation 
direction  z.  The  geometry  of  the  problem  is  illustrated  in  Figure  1. 
For  reasons  of  computational  efficiency  to  be  demonstrated  later,  (2) 


is  in  some  cases  preferable  to  a more  standard  approach  > in  which 

ik.z 

the  phase  factor  is  taken  to  be  e 


0 . Substitution  of  (2)  into  (l) 


gives 


where 


V2u  + 2ik(z)  ?)  u + u(i3  k + k * x ) = 0 
Z z O 1 


(M 


X1  (£ ) = x(£)  - x0'z)  • 

One  further  transformation  is  useful  before  introducing  the  parabolic 


approximation:  let 


f (r ) = u(r)  (k/k  )2 


(6) 


u(r  ) (1  + xQ'z)' 


. i. 


Thus 


E (r ) = f (r ) (1  + x„  ( z ) ) 4 e 


J k(z  )dz 


(7) 


Then  (k)  becomes 

^i2f  + k^  g 2 (fk~2)  + 2ika_f  + fk  2X  - 0 , (8) 

z z o i 

Where  i refers  to  the  two  Cartesian  coordinates  x and  y orthogonal  to 
z.  It  is  important  to  note  that  no  approximations  have  been  introduced 
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between  (l)  and  (3). 

2.  THE  PARABOLIC  APPROXIMATION 


If  we  assume  that  f is  substantially  constant  over  a distance  of 
one  wavelength 


X * 2n/ko  , 

(9) 

and  that  !y  1 « 1 , 

1 

(10) 

the  second  z derivative  in  (8)  becomes  negligible  in  comparison  with 
the  first  z derivative.  To  define  what  "negligible"  means,  let  us 
assume  that  f varies  on  a length  scale  L » X and  assign  constant 
values  to  k and  y • We  Fourier  transform  f according  to 


The  transform  of  (8)  then  gives 

•P  2 • P 2 • 2k  P + k 2V  * 0 . (12) 

1 Z Z O *•],  v 

Because  of  the  slow  variation  assumed  for  f,  we  can  assign  a limiting 

value  to  P : 
j. 

P^2  < 1/L2  . (13) 

If  we  take  k k , (12)  and  (13 ) give  a limiting  value  for  P : 
o z 


(lM 


The  right  side  of  (lh- ) is  small  by  assumption,  so 

Since  the  quadratic  term  can  be  dropped  from  (lM,  the  second  z 


(15) 
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derivative  may  be  dropped  from  (8): 


r 


a f = — (v  2f  + k 2x  f ) . 

Z £lK  ± O 1 

This  is  the  "parabolic  equation"  for  f.  The  fractional  error  in 
resulting  from  the  parabolic  approximation  is  on  the  order  of 


(16) 

f 


3.  NUMERICAL  INTEGRATION  OF  THE  PARABOLIC  EQUATION 


Note  that  (1 6)  represents  an  initial  value  problem  with  regard 
to  z integration.  We  can  specify  f at  some  initial  z to  describe  an 
incident  wave,  and  use  (l6)  to  step  through  a region  of  arbitrarily 
varying  x • Previous  numerical  application  of  the  parabolic  equation 
to  the  ionospheric  scintillation  problem  has  been  primarily  oriented 
toward  propagation  into  free  space  of  a wave  which  has  just  emerged 
from  a random  phase  screen4.  In  this  application,  one  Fourier  trans- 
forms the  initial  f(x,y,o),  advances  the  phase  of  each  spectral 
component  in  an  appropriate  way,  and  then  transforms  back  to  get 
f(x,y,z).  However,  by  integrating  (16 ) within  the  scattering  medium, 
one  can  follow  the  growth  of  amplitude  and  phase  variations  in  an 
arbitrarily  thick  slab. 


A.  FINITE  DIFFERENCE  SCHEME 

We  propose  a leapfrog  scheme  to  integrate  (l6),  with  real  and 
imaginary  parts  of  f defined  on  alternate  z levels  (see  Figure  2). 

Let  us  express  f explicitly  in  terms  of  real  and  imaginary  parts: 


f(r)  - fr'r)  + ifi(r)  (17) 

where  f and  f ^ are  real.  Then  (1 6)  gives  the  coupled  linear  system: 

S f * zf  h 2f,  + k 2y  f,') 

zr2klxio*1iJ  (lo ) 

3 f,  + k2vf\ 

zt  2k  y x r o^iry- 

This  system  may  be  integrated  on  a staggered  mesh,  with  fr  and 
defined  on  alternate  z levels.  The  three  functions  fr,  f^,  and  ^ 
are  discretized  in  (x,y,z)  space  with  integer  indices  (I,J,K).  The 
mesh  spacings  in  the  x and  y directions  are  the  same  value,  fx;  the 
spacing  between  the  z planes  is  6z . We  shall  use  the  following  second 
order  difference  formulas: 


7^2  f (I, J ,K ) [f(I  + 1, J,K ) + f(I-l,J,K) 

+ f(I,J  + 1 ,K ) + f (I , J-l  ,K ) - 4f(I,J,K)]/«X2 


(19) 


B.  STABILITY  ANALYSIS 

If  an  additive  error  € = e + ie  ^ is  introduced  into  a discre- 
tized solution  f,  then  e itself  must  satisfy  (l8).  This  error  may  be 
Fourier  analyzed,  and  stability  analysis  applied  to  each  spectral 
component.  We  take  ^ constant,  set  k = kQ  for  simplicity,  and  take 
i(PX  + PY  + PZ) 

e(r)  <*  e X y 2 . (20 ) 

For  a specific  choice  of  (Px,P  ,P2),  the  derivative  operators  Sz 

and  7 2 may  be  treated  as  algebraic  multipliers.  Substitution  of  (20 ) 
x 

into  (1 6)  gives 


5 


(21) 


where 


) 


and 


iP  5 -iP  *z 
z z z 

e - e 


2'iz 

is  in  (Pz*z)/fz  , 

2 C2  _ cos  'px5x) 


COS 'Py*X)] 


(22) 


(25) 


For  investigating  stability  of  integrating  forward  in  the  z direction, 
we  set 


iP  *z 

V = e Z (2k) 

The  stability  criterion  will  be  |yj  ■s.  1.  Now  (2l)  gives  the  quadratic 
equation 

V - Y'1  = iu  , (25) 

where 


The  solution  to  (25)  is 

Y = i|  ± 7l  - u A (27  ) 

2 

When  u A > 1,  one  of  the  roots  gives  |y|  >1,  and  the  scheme  is 
2 

unstable.  When  M-  A s 1,  the  radical  in  (27)  is  real,  and 

Ivl2  = + (1  - ,2A)  = 1.  (28) 

Thus  the  stability  criterion  is 

In!  £ 2 (29) 

The  greatest  possible  value  of  |y,|  occurs  when  7^2  = -8/6X2,  and 
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X ^ < 0.  Thus  all  modes  are  stable  when  the  following  is  satisfied: 


STABILITY 

REQUIREMENT: 


2k 


5z  £ 


8/5x^  + k^lx  I 


(30) 


This  criterion  demonstrates  the  advantage  of  separating  the  background 

ik<}Z 

X from  the  fluctuations  x . If  we  had  used  e as  a phase  factor 


in  (2),  then  x in  (30)  would  be  replaced  by  xQ  + 


In  typical 


applications,  this  would  result  in  a smaller  stepsize  sz. 

In  addition  to  the  stability  criterion  (50)  we  wish  also  to 
impose  a resolution  requirement  that  a change  by  only  a "small  amount" 
from  one  z level  to  the  next.  By  this  we  mean  that  the  function  f in 
(lo)  should  undergo  a phase  change  less  than  one  radian  over  a distance 
26z: 


|2pz5z(  < 1 , 

where  the  limiting  form  (15 ) is  used  for  p 


(31) 


Thus  we  have 


RESOLUTION 

REQUIREMENT: 


Az  < k 


-1  (w + 


(32) 


We  have  made  changes  in  going  to  (15)  to  (55 ) to  account  for  the  fact 
that  x may  be  positive  or  negative.  The  ratio  of  the  largest  6z  per- 
mitted by  the  resolution  condition  to  the  largest  value  allowed  by  the 
stability  condition  is 
AZ. 

\ u ..  /I  , 

= r-  ! 

(33) 


AZ, 


8/ax2  + kQ2  |Xi 
2 (l/L2  + k 


'V 


where  5zn  and  az  are  taken  from  (52)  and  (51)  respectively.  For  small 
R ^ 


* 


X , Chis  ratio  tends  to  4L2/ax2.  This  value  is  by  assumption  greater 
than  one  for  all  cases  that  can  be  considered  within  the  framework  of 
a discrete  point  representation.  When  the  term  dominates,  (33) 


tends  toward  the  value  -/2.  To  insure  stability  and  adequate  resolu- 
tion, we  take  the  lesser  of  (30  ) and  32 ) for  hz . This  can  be 
summarized  as  follows: 

When  kQ2  | s 8/^ix2  - 2/L^  use  (30  ) 

" > " use  (32). 

We  have  treated  I as  a constant  in  writing  (JO)  - (33).  In 

practice,  one  would  use  the  maximum  value  at  a given  z level  to  deter- 
mine the  stepsize  sz . 

C.  BOUNDARY  CONDITIONS 

In  the  examples  to  be  presented  in  the  next  section  we  have  used 
a very  simple  (although  arbitrary)  algorithm  for  defining  the 
Laplacian  operator  7 2 at  the  edges  of  the  mesh.  Namely,  we  have  set 
the  normal  second  derivative  equal  to  zero  at  all  locations  one  cell 
from  the  boundary.  We  confine  the  calculation  to  the  interior  mesh. 
Suppose  that  the  (x,y)  mesh  consists  of  NX  by  NY  points.  Then  for  all 
J and  K we  impose  the  condition 

f (1, J,K)  = 2f(2,J,K)  - f(3,J,K)  (3*0 

f (NX, J ,K ) = 2f (NX-1, J,K ) - f (NX-2 , J ,K ) . 

A straightforward  extension  of  (3*0  defines  f for  J = 1 and  NY.  In 
practice,  we  use  (l8 ) and  (19 ) to  update  f on  the  interior  mesh 
I = 2 to  NX-1,  J - 2 to  NY-1,  and  define  f on  the  boundaries  with  (34). 
We  have  so  far  not  noticed  any  stability  or  accuracy  problems  resulting 

8 


1 


from  (34 ) . 


D.  RESULTS  FROM  THE  PEP  CODE 


We  have  two  comparisons  to  present  between  numerical  solutions  to 
(18)  and  published  results  obtained  by  other  means.  We  have  carried 
out  the  numerical  procedure  defined  in  (18 ) and  (19)  and  included  cer- 
tain diagnostic  features  in  a "Parabolic  Equation  Propagation"  code, 
whose  acronym  is  conveniently  PEP.  Values  for  the  real  and  imaginary 
parts  of  f are  assigned  at  some  initial  z level.  The  PEP  code  steps 
(l8)  forward,  and  at  specified  intervals  computes  "scintillation 
indices,"  phase  variances  in  radians,  and  power  variances  in  decibels. 
The  scintillation  index  of  greatest  interest  is 

. V<m4>  - 

V*  ’ (55) 

where  the  bracket  refers  to  an  average  over  x and  y,  with  z fixed. 
The  phase  and  power  variances  re  respectively 


= <$2>  - <$> 


where 

$ = Tan-1  J and 

a 2(z)  = 100  < (Log  lf!2)2>  - 100  < Log  |f|2>  . (37) 

X 10  10 

The  notation  conflict  between  the  power  variance  a and  the  suscepti- 

X 

bility  i(  is  a result  of  established  conventions  in  ionospheric  and 
plasma  physics  literature. 


Focusing  by  a Gaussian  Lens 


First  we  consider  a wave  incident  on  a two  dimensional  lens 
x,  z)  whose  thickness  in  terms  of  phase  delay  is 

$(x,o)  = 10  exp  ( - (x/j>J2)  (58) 

radians.  This  problem  has  been  examined  numerically  by  Buckley5  using 
a small  angle  diffraction  approach.  Our  approach  is  to  solve  the 
initial  value  problem  (l6)  with  y =0  everywhere  below  the  lens 
(z  < 0),  and  with  initial  conditions 

f(x,o)  - e ’ . (39) 

The  two  dimensional  problem  was  integrated  according  to  the  prescrip- 
tion (l8)  - (19)  with  the  J index  suppressed.  Numerical  values  used 
for  the  integration  were:  fix  = .05,  fiz  = 0.31^159,  X = .01,  l = 1, 

NZ  = 102,  and  NX  = 82 . Results  for  the  "power"  |E(x,z)|2  = |f(x,z)|2are 
in  Figure  J . The  qualitative  agreement  with  Buckley's  result5  is 
quite  satisfactory.  Z in  Figure  3 is  a dimensionless  coordinate  which 
can  be  used  to  scale  equation  (l6): 

Z = z /k£2  ( lO  ) 


2.  Diffraction  from  a Random  Phase  Screen 


Let  us  consider  a three  dimensional  problem  in  which  the  initial 
phase  $(x,y,o)  is  a random  function  of  x and  y with  Gaussian  auto- 
correlation : 


-(?2  + lf)/r02 

<$(x,y,o)  $(x  + e,y  + T1,0)>  = $Q2e  • (^l) 


Such  a function  can  be  generated  numerically  by  Fourier  transforming 
the  desired  spatial  autocorrelation  function,  multiplying  each  Fourier 


1 
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i,- 

component  by  a random  phase  factor  , and  then  transforming  back 

into  xy  space.  In  this  prescription,  o is  an  arbitrarily  generated 

number  ranging  from  0 to  2n.  To  insure  that  $ is  real,  one  can  assign 

cd  arbitrarily  in  one  half  of  the  transform  space  (k  , k ) and  then  set 

x'  y 

»(-k  , -k  ) ■ - «(k  , k ) . 
x'  y x'  y 

We  have  generated  a function  satisfying  (4l)  to  use  as  an  initial 
condition  for  (l6).  We  initialize  the  calculation  with 

f(x,y,o)  = g1*  'X>y^°  ) (42) 

and  integrate  forward  into  free  space  (x  * 0).  Four  curves  showing 
computed  values  of  S versus  propagation  distance  for  different  values 

4 

of  appear  in  Figure  4.  The  dots  are  from  the  calculation,  and  the 
solid  curves  are  the  results  of  Mercier6  plotted  in  the  format  of 
Briggs  and  Parkin7.  The  results  are  in  good  agreement  with  the  solid 
curves.  The  data  used  to  generate  these  results  were  ^x  = .002, 

5z  * .015,  X = .0001,  rQ  = .01,  NX  = 27,  NY  = 27,  and  NZ  = 102. 

SUMMARY 


We  have  described  a set  of  coupled  linear  equations  (l8  ) which 
result  from  the  wave  equation  (l ) in  the  parabolic  approximation. 
Equations  (l8)  differ  from  the  usual  "parabolic  equation"  only  in 
that  the  index  of  refraction  is  allowed  to  vary  slightly  along  an 
assumed  propagation  direction  as  part  of  the  background  or  unscattered 
state.  A rough  estimate  of  the  error  introduced  by  the  parabolic 
approximation  is  given  after  (l6).  We  have  constructed  the  PEP  code 
for  finite  difference  solution  of  (l8 ) using  the  centered  leapfrog 
scheme  (19)*  The  stability  criterion  for  the  scheme  is  (30). 
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Results  from  the  PEP  code  are  in  good  agreement  with  those 
obtained  by  different  means  for  refraction  due  to  a Gaussian  lens  and 
for  refraction  due  to  a random  phase  screen.  These  two  examples 
involve  propagation  into  a vacuum  of  an  initially  specified  wave. 
Results  for  propagation  within  a scattering  medium  will  be  given  else- 
where . 
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INCIDENT  WAVE 


SCATTERED  WAVE 

SCATTERING  VOLUME 


Fig.  1 — A wave  incident  along  the  z direction  is  scattered  by 
inhomogeneities  in  the  scattering  volume 


A B 


Fig.  2 — Leapfrog  difference  scheme.  Step  A updates  the  real  part  of  f 
from  plane  1 to  plane  3 according  to  (18)  and  (19).  Then  Step  B up- 
dates the  imaginary  part  of  f from  plane  2 to  plane  4. 
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A 


0- 


I A 2 / r02 


Fig.  4 — Scintillation  index  as  a function  of  distance  from  a random  phase  screen. 
The  curves,  labelled  with  their  initial  phase  variances  in  radians,  are  from  Briggs  and 
Parkin7 . The  dots  are  results  from  the  PIP  code. 
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